Report with geocomputation examples for point, line and polygon data

This document provides the methodlogy and the main results regarding the spatialization of pollutation inventory.

# Function for preparing point data for spatialization 
# Parameters:
#  - source.list - list with source data and total inventory
#  - distribute - TRUE or FALSE, depend on that if you want to distribute total inventory emissions to sources, when difference exists
corsum2sf <- function(source.list, distribute = FALSE){
  source.list$sources$points[, vars] <- source.list$sources$points[, vars] %>% dplyr::mutate_all(~replace(., is.na(.), 0))
  source.list[[2]][[2]][, vars] <- source.list[[2]][[2]][, vars] %>% dplyr::mutate_all(~replace(., is.na(.), 0))
  
  points.sum <- source.list$sources$points %>% 
    dplyr::select(., vars) %>% 
    apply(., 2, sum) %>% 
    t(.) %>% 
    as.data.frame() %>%
    dplyr::mutate_if(is.numeric, round, 2)
  
  points.total <- source.list[[2]][[2]][, vars] %>% 
    dplyr::mutate_if(is.numeric, round, 2) %>%
    as.data.frame()
  
  if(!identical(points.sum, points.total) & distribute == TRUE){
    d <- (points.total - points.sum)[1, ]
    zero.ind <- source.list$sources$points[, vars] == 0
    w <- replace(source.list$sources$points[, vars], zero.ind, 1) %>% 
      apply(., 2, function(x) x/sum(x)) %>%
      as.data.frame() %>%
      dplyr::mutate_all(~replace(., is.na(.), 0))
    
    cor.data <- as.matrix(w) %*% diag(d) + source.list$sources$points[, vars]
    source.list$sources$points[, vars] <- cor.data
    source.sf <- st_as_sf(x = source.list$sources$points, coords = c("Longitude", "Latitude"),  crs = "+proj=longlat +datum=WGS84")
  }else{
    source.sf <- st_as_sf(x = source.list$sources$points, coords = c("Longitude", "Latitude"),  crs = "+proj=longlat +datum=WGS84")
  }
  return(source.sf)
}

# Function for preparing polygon data for spatialization 
# Parameters:
#  - source.list - list with source data and total inventory
#  - distribute - TRUE or FALSE, depend on that if you want to distribute total inventory emissions to sources, when difference exists

corsum2sf_polygon <- function(source.list, distribute = FALSE){
  source.list$sources$polygon[, vars] <- source.list$sources$polygon[, vars] %>% dplyr::mutate_all(~replace(., is.na(.), 0)) %>% st_drop_geometry()
  source.list[[2]][[2]][, vars] <- source.list[[2]][[2]][, vars] %>% dplyr::mutate_all(~replace(., is.na(.), 0))
  
  polygon.sum <- source.list$sources$polygon %>% 
    st_drop_geometry() %>%
    dplyr::select(., vars) %>% 
    apply(., 2, sum) %>% 
    t(.) %>% 
    as.data.frame() %>%
    dplyr::mutate_if(is.numeric, round, 2)
  
  polygon.total <- source.list[[2]][[2]][, vars] %>% 
    dplyr::mutate_if(is.numeric, round, 2) %>%
    as.data.frame()
  
  if(!identical(polygon.sum, polygon.total) & distribute == TRUE){
    d <- (polygon.total - polygon.sum)[1, ]
    zero.ind <- source.list$sources$polygon[, vars] %>% st_drop_geometry() == 0
    w <- replace(source.list$sources$polygon[, vars] %>% st_drop_geometry(), zero.ind, 1) %>% 
      apply(., 2, function(x) x/sum(x)) %>%
      as.data.frame() %>%
      dplyr::mutate_all(~replace(., is.na(.), 0))
    
    cor.data <- as.matrix(w) %*% diag(d) + source.list$sources$polygon[, vars]%>% st_drop_geometry()
    source.list$sources$polygon[, vars] <- cor.data
    source.sf <- source.list$sources$polygon
  }else{
    source.sf <- source.list$sources$polygon
  }
  return(source.sf)
}

# Function for preparing line data for spatialization 
# Parameters:
#  - source.list - list with source data and total inventory
#  - distribute - TRUE or FALSE, depend on that if you want to distribute total inventory emissions to sources, when difference exists

corsum2sf_lines <- function(source.list, distribute = FALSE){
  source.list$sources$lines[, vars] <- source.list$sources$lines[, vars] %>% dplyr::mutate_all(~replace(., is.na(.), 0)) %>% st_drop_geometry()
  source.list[[2]][[2]][, vars] <- source.list[[2]][[2]][, vars] %>% dplyr::mutate_all(~replace(., is.na(.), 0))
  
  lines.sum <- source.list$sources$lines %>% 
    st_drop_geometry() %>%
    dplyr::select(., vars) %>% 
    apply(., 2, sum) %>% 
    t(.) %>% 
    as.data.frame() %>%
    dplyr::mutate_if(is.numeric, round, 2)
  
  lines.total <- source.list[[2]][[2]][, vars] %>% 
    dplyr::mutate_if(is.numeric, round, 2) %>%
    as.data.frame()
  
  if(!identical(lines.sum, lines.total) & distribute == TRUE){
    d <- (lines.total - lines.sum)[1, ]
    zero.ind <- source.list$sources$lines[, vars] %>% st_drop_geometry() == 0
    w <- replace(source.list$sources$lines[, vars] %>% st_drop_geometry(), zero.ind, 1) %>% 
      apply(., 2, function(x) x/sum(x)) %>%
      as.data.frame() %>%
      dplyr::mutate_all(~replace(., is.na(.), 0))
    
    cor.data <- as.matrix(w) %*% diag(d) + source.list$sources$lines[, vars]%>% st_drop_geometry()
    source.list$sources$lines[, vars] <- cor.data
    source.sf <- source.list$sources$lines
  }else{
    source.sf <- source.list$sources$lines
  }
  return(source.sf)
}

# Function for spatial data visualisation at web maps
# Parameters:
#    1. sf.sources -  sf object with sources to be spatialised
#    2. sf.spatialised - sf object with polygons after spatialisation
spatialised.mapview <- function(sf.sources, layer.name.1 = "", sf.spatialised, layer.name.2 = "", vars = vars, source.lines = FALSE){
  sf.spatialised$Spatialised <- NA
  sf.spatialised[ ,vars] %<>% st_drop_geometry() %>% dplyr::mutate_all(.,as.double)
  sf.spatialised$Spatialised[sf.spatialised$NOx == 0 & sf.spatialised$SO2 == 0 & sf.spatialised$PM10 == 0 & sf.spatialised$PM2.5 == 0 & sf.spatialised$NMVOC == 0 & sf.spatialised$NH3 == 0] <- 0
  sf.spatialised$Spatialised[sf.spatialised$NOx !=0 | sf.spatialised$SO2 !=0 | sf.spatialised$PM10 !=0 | sf.spatialised$PM2.5 !=0 | sf.spatialised$NMVOC !=0 | sf.spatialised$NH3 != 0] <- 1
  
  if(source.lines == FALSE){
    web_map <- mapview(sf.spatialised, layer.name = layer.name.2, zcol = "Spatialised") + mapview(sf.sources, layer.name = layer.name.1, col.regions = "red")
  }else {
    web_map <- mapview(sf.sources, layer.name = layer.name.1, color = "red") +  mapview(sf.spatialised, layer.name = layer.name.2, zcol = "Spatialised")
  }
  return(web_map)
}

Example - point data

# Reading source file and grid
source.file = "Pollutant inventory spatialized-d30102019.xlsx"
source.sheet =  "1A1-Energy"
header <- readxl::read_xlsx(path = source.file, range = "D8:S8", sheet = source.sheet) %>% names()

## New names:
## * `` -> ...10

vars <- header[1:6]
grid.5km <- readOGR("Grid/Polygons_5km_UTM_34N.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\R_projects\Spatialization\Grid\Polygons_5km_UTM_34N.shp", layer: "Polygons_5km_UTM_34N"
## with 3124 features
## It has 1 fields
## Integer64 fields read as strings:  ID

sf.grid.5km <- st_as_sf(grid.5km)

# Reading source data
source.1A1a <- list(sources = list(points = NA, lines = NA, polygon = NA), total = list(spatialize = NA, inventory = NA))

source.1A1a$sources$points <- readxl::read_xlsx(path = source.file, range = "D9:S37", sheet = source.sheet, col_names = header)
source.1A1a$total$spatialize <- readxl::read_xlsx(path = source.file, range = "D40:I40", sheet = source.sheet, col_names = vars)
source.1A1a$total$inventory <- readxl::read_xlsx(path = source.file, range = "D46:I46", sheet = source.sheet, col_names = vars)

# Preparing data for final spatialization
sf.1A1a <- corsum2sf(source.1A1a) %>%
  st_transform(crs = "+init=epsg:32634")

sf.1A1a %>% 
  st_drop_geometry() %>% 
  dplyr::mutate_if(is.numeric, round, 2) %>%
  datatable(., caption = 'Table 1: sf.1A1a',
            rownames = FALSE, escape = FALSE, selection = "single",
            extensions = c('Buttons'),
            class = 'white-space: nowrap',
            options = list(
              pageLength = 5,
              dom = 'Bfrtip',
              buttons = list('pageLength'),
              searchHighlight = TRUE,
              scrollX = TRUE,
              scrollY = TRUE
            ))

sum.1A1a <- sf.1A1a %>% 
  st_drop_geometry() %>%
  dplyr::select(., vars) %>% 
  apply(., 2, sum) %>% 
  t(.) %>% 
  as.data.frame() %>%
  dplyr::mutate_if(is.numeric, round, 2)
total.1A1a <- source.1A1a[[2]][[2]][, vars] %>% 
  mutate_all(~replace(., is.na(.), 0)) %>%
  dplyr::mutate_if(is.numeric, round, 2) %>%
  as.data.frame()

data.frame(sum = c("spatialize", "total", "diff"), rbind(sum.1A1a, total.1A1a, data.frame(sum.1A1a == total.1A1a)-1)) %>%
  datatable(., caption = 'Table 2: Summary differences',
            options = list(pageLength = 5)
  )

p.1A1a <- sf.grid.5km %>% # Spatialization by grid cells
  st_join(sf.1A1a) %>%
  group_by(ID) %>%
  summarize(NOx = sum(NOx, na.rm = TRUE), # Sum-up by grid cells
            SO2 = sum(SO2, na.rm = TRUE),
            PM10 = sum(PM10, na.rm = TRUE),
            PM2.5 = sum(PM2.5, na.rm = TRUE),
            NMVOC = sum(NMVOC, na.rm = TRUE),
            NH3 = sum(NH3, na.rm = TRUE)) %>% 
  mutate(ID = as.numeric(ID))

sum.p.1A1a <- p.1A1a %>% # Computational control
  st_drop_geometry() %>%
  dplyr::select(., vars) %>% 
  apply(., 2, sum) %>% 
  t(.) %>% 
  as.data.frame() %>%
  dplyr::mutate_if(is.numeric, round, 2)
data.frame(sum = c("spatialized", "total", "diff"), rbind(sum.p.1A1a, total.1A1a, data.frame(sum.p.1A1a == total.1A1a)-1)) %>%
  datatable(., caption = 'Table 3: Summary differences after spatialisation',
            options = list(pageLength = 5)
  )

Example - line data

# Reading source file and grid
source.file = "Pollutant inventory spatialized-d30102019.xlsx"
source.sheet =  "1A3-Transport"
header <- readxl::read_xlsx(path = source.file, range = "D9:S9", sheet = source.sheet) %>% names()
vars <- header[1:6]
grid.5km <- readOGR("Grid/Polygons_5km_UTM_34N.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\R_projects\Spatialization\Grid\Polygons_5km_UTM_34N.shp", layer: "Polygons_5km_UTM_34N"
## with 3124 features
## It has 1 fields
## Integer64 fields read as strings:  ID

sf.grid.5km <- st_as_sf(grid.5km) 

# Reading source data
source.1A3c <- list(sources = list(points = NA, lines = NA, polygon = NA), total = list(spatialize = NA, inventory = NA))

source.1A3c$total$spatialize <- readxl::read_xlsx(path = source.file, range = "D159:I159", sheet = source.sheet, col_names = vars)
source.1A3c$total$inventory <- readxl::read_xlsx(path = source.file, range = "D160:I160", sheet = source.sheet, col_names = vars)

# Reading data
railways <- readOGR("Data/pruge/Pruge_osm_32634.shp", 
                    use_iconv=TRUE,  
                    encoding = "UTF-8",
                    stringsAsFactors = FALSE)

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\R_projects\Spatialization\Data\pruge\Pruge_osm_32634.shp", layer: "Pruge_osm_32634"
## with 4921 features
## It has 4 fields

# Transformation to sf:: class and transformation to UTM reference coordinate system
sf_railways <- st_as_sf(railways) %>%
  st_transform(crs = "+init=epsg:32634") 

# Make new columns
sf_railways[,vars] <- NA
# Intersection with grid cells
sf_railways.int <- st_intersection(sf_railways, sf.grid.5km) %>%
  select(.,vars, osm_id) %>%
  mutate(Length = st_length(.))

# Preparing data for final spatialization
source.1A3c$sources$lines <- sf_railways.int
sf.1A3c <- corsum2sf_lines(source.1A3c, distribute = FALSE)

sf.1A3c %>% 
   st_drop_geometry() %>% 
   dplyr::mutate_if(is.numeric, round, 2) %>%
   datatable(., caption = 'Table 4: sf.1A3c',
             rownames = FALSE, escape = FALSE, selection = "single",
             extensions = c('Buttons'),
             class = 'white-space: nowrap',
             options = list(
               pageLength = 5,
               dom = 'Bfrtip',
               buttons = list('pageLength'),
               searchHighlight = TRUE,
               scrollX = TRUE,
               scrollY = TRUE
             ))

sum.1A3c <- sf.1A3c %>% 
  st_drop_geometry() %>%
  dplyr::select(., vars) %>% 
  apply(., 2, sum) %>% 
  t(.) %>% 
  as.data.frame() %>%
  dplyr::mutate_if(is.numeric, round, 2)
total.1A3c <- source.1A3c[[2]][[2]][, vars] %>% 
  mutate_all(~replace(., is.na(.), 0)) %>%
  dplyr::mutate_if(is.numeric, round, 2) %>%
  as.data.frame()

data.frame(sum = c("spatialize", "total", "diff"), rbind(sum.1A3c, total.1A3c, data.frame(total.1A3c - sum.1A3c))) %>%
  datatable(., caption = 'Table 5: Summary differences',
            options = list(pageLength = 5)
  )

sum_Length <- sum(sf.1A3c$Length)
diff.1A3c <- data.frame(total.1A3c - sum.1A3c)
sf.1A3c <- sf.1A3c %>% # Calculate weights and distribute values
  mutate(NOx = ((diff.1A3c$NOx/sum_Length)*Length),
         SO2 = ((diff.1A3c$SO2/sum_Length)*Length),
         PM10 = ((diff.1A3c$PM10/sum_Length)*Length),
         PM2.5 = ((diff.1A3c$PM2.5/sum_Length)*Length),
         NMVOC = ((diff.1A3c$NMVOC/sum_Length)*Length),
         NH3 = ((diff.1A3c$NH3/sum_Length)*Length))
sf.1A3c %<>% select(vars)

p.1A3c <- sf.grid.5km %>% # Spatialization by grid cells
  st_join(sf.1A3c, join = st_contains) %>% 
  group_by(ID) %>%
  summarize(NOx = sum(NOx, na.rm = TRUE), # Sum-up by grid cells
            SO2 = sum(SO2, na.rm = TRUE),
            PM10 = sum(PM10, na.rm = TRUE),
            PM2.5 = sum(PM2.5, na.rm = TRUE),
            NMVOC = sum(NMVOC, na.rm = TRUE),
            NH3 = sum(NH3, na.rm = TRUE)) %>% 
  mutate(ID = as.numeric(ID))

sum.p.1A3c <- p.1A3c %>% # Computational control
  st_drop_geometry() %>%
  dplyr::select(., vars) %>% 
  apply(., 2, sum) %>% 
  t(.) %>% 
  as.data.frame() %>%
  dplyr::mutate_if(is.numeric, round, 2)
data.frame(sum = c("spatialized", "total", "diff"), rbind(sum.p.1A3c, total.1A3c, data.frame(sum.p.1A3c == total.1A3c)-1)) %>%
  datatable(., caption = 'Table 6: Summary differences after spatialisation',
            options = list(pageLength = 5)
  )

# Export data to geopackage file
# st_write(p.1A3c, dsn="Products/1A3 - Transport/1A3c.gpkg", layer='1A3c')

Example - polygon data

# Reading source file and grid
source.file = "Pollutant inventory spatialized-d30102019.xlsx"
source.sheet =  "1A4-Residential-Tertiary"
header <- readxl::read_xlsx(path = source.file, range = "D8:S8", sheet = source.sheet) %>% names()
vars <- header[1:6]
grid.5km <- readOGR("Grid/Polygons_5km_UTM_34N.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\R_projects\Spatialization\Grid\Polygons_5km_UTM_34N.shp", layer: "Polygons_5km_UTM_34N"
## with 3124 features
## It has 1 fields
## Integer64 fields read as strings:  ID

sf.grid.5km <- st_as_sf(grid.5km) 

# Reading source data
source.1A4ai <- list(sources = list(points = NA, lines = NA, polygon = NA), total = list(spatialize = NA, inventory = NA))

source.1A4ai$total$spatialize <- readxl::read_xlsx(path = source.file, range = "D18:I18", sheet = source.sheet, col_names = vars)
source.1A4ai$total$inventory <- readxl::read_xlsx(path = source.file, range = "D22:I22", sheet = source.sheet, col_names = vars)

# Reading data - CLC - Corine Land Cover
clc_18 <- readOGR("Data/clc/CLC18_RS.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\R_projects\Spatialization\Data\clc\CLC18_RS.shp", layer: "CLC18_RS"
## with 34794 features
## It has 6 fields

sf_clc18 <- st_as_sf(clc_18)
# Subset data to get requested CLC class
sf_clc18_urb <- subset(sf_clc18, CODE_18 == "111" | CODE_18 == "112") %>% # CLC urban zones
  st_transform(crs = "+init=epsg:32634")

# Readinf data - munucipalites polygons
Sys.setlocale(locale = 'Serbian (Latin)')

## [1] "LC_COLLATE=Serbian (Latin)_Bosnia and Herzegovina.1250;LC_CTYPE=Serbian (Latin)_Bosnia and Herzegovina.1250;LC_MONETARY=Serbian (Latin)_Bosnia and Herzegovina.1250;LC_NUMERIC=C;LC_TIME=Serbian (Latin)_Bosnia and Herzegovina.1250"

opstine <- readOGR("Data/opstine/gadm36_SRB_2.shp", 
                   use_iconv=TRUE,  
                   encoding = "UTF-8")

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\R_projects\Spatialization\Data\opstine\gadm36_SRB_2.shp", layer: "gadm36_SRB_2"
## with 161 features
## It has 13 fields

sf_opstine <- st_as_sf(opstine)

# Read data from  Official annual report from the public company “Toplane Srbije” - Total_heating values
toplane <- readxl::read_xls(path = "Data/toplane/Toplane_2015.xls") %>%
  mutate(GRAD = str_to_title(GRAD))
sf_opstine$Toplane <- toplane$`Ukupna grejna povrsina (m2)`[match(sf_opstine$NAME_2, toplane$GRAD)] 
sf_opstine %<>% 
  st_transform(crs = "+init=epsg:32634") %>% 
  mutate_all(~replace(., is.na(.), 0))

# Join data - from munucipalities ('opstine') to CLC urban class
sf_clc18_urb <- st_join(sf_clc18_urb, sf_opstine, largest = TRUE) 
sf_clc18_urb %<>% dplyr::select(.,Toplane, NAME_2)
sf_clc18_urb[,vars] <- NA

# Intersection with grid cells
sf_clc18_urb.int <- st_intersection(sf_clc18_urb, sf.grid.5km) %>% 
  filter(!is.na(Toplane))

# Preparing data for final spatialization
source.1A4ai$sources$polygon <- sf_clc18_urb.int
sf.1A4ai <- corsum2sf_polygon(source.1A4ai, distribute = FALSE) %>%
  st_transform(crs = "+init=epsg:32634")

sf.1A4ai %>% 
  st_drop_geometry() %>%
  dplyr::rename(Total_heating_area = Toplane, Manucipality = NAME_2) %>%
  dplyr::mutate_if(is.numeric, round, 2) %>%
  datatable(., caption = 'Table 1: sf.1A4ai',
            rownames = FALSE, escape = FALSE, selection = "single",
            extensions = c('Buttons'),
            class = 'white-space: nowrap',
            options = list(
              pageLength = 5,
              dom = 'Bfrtip',
              buttons = list('pageLength'),
              searchHighlight = TRUE,
              scrollX = TRUE,
              scrollY = TRUE
            ))

sum.1A4ai <- sf.1A4ai %>% 
  st_drop_geometry() %>%
  dplyr::select(., vars) %>% 
  apply(., 2, sum) %>% 
  t(.) %>% 
  as.data.frame() %>%
  dplyr::mutate_if(is.numeric, round, 2)
total.1A4ai <- source.1A4ai[[2]][[2]][, vars] %>% 
  mutate_all(~replace(., is.na(.), 0)) %>%
  dplyr::mutate_if(is.numeric, round, 2) %>%
  as.data.frame()

data.frame(sum = c("spatialize", "total", "diff"), rbind(sum.1A4ai, total.1A4ai, data.frame(total.1A4ai - sum.1A4ai))) %>%
  datatable(., caption = 'Table 2: Summary differences',
            options = list(pageLength = 5)
  )

sum_s <- sum(sf.1A4ai$Toplane)
diff.1A4ai <- data.frame(total.1A4ai - sum.1A4ai)
sf.1A4ai <- sf.1A4ai %>% # Calculate weights and distribute values
  mutate(NOx = ((diff.1A4ai$NOx/sum_s)*Toplane),
         SO2 = ((diff.1A4ai$SO2/sum_s)*Toplane),
         PM10 = ((diff.1A4ai$PM10/sum_s)*Toplane),
         PM2.5 = ((diff.1A4ai$PM2.5/sum_s)*Toplane),
         NMVOC = ((diff.1A4ai$NMVOC/sum_s)*Toplane),
         NH3 = ((diff.1A4ai$NH3/sum_s)*Toplane))
sf.1A4ai %<>% select(vars)

p.1A4ai <- sf.grid.5km %>% # Spatialization by grid cells
  st_join(sf.1A4ai, join = st_contains) %>% 
  group_by(ID) %>%
  summarize(NOx = sum(NOx, na.rm = TRUE), # Sum-up by grid cells
            SO2 = sum(SO2, na.rm = TRUE),
            PM10 = sum(PM10, na.rm = TRUE),
            PM2.5 = sum(PM2.5, na.rm = TRUE),
            NMVOC = sum(NMVOC, na.rm = TRUE),
            NH3 = sum(NH3, na.rm = TRUE)) %>% 
  mutate(ID = as.numeric(ID))

sum.p.1A4ai <- p.1A4ai %>% # Computational control
  st_drop_geometry() %>% 
  dplyr::select(., vars) %>% 
  apply(., 2, sum) %>% 
  t(.) %>% 
  as.data.frame() %>%
  dplyr::mutate_if(is.numeric, round, 2)
data.frame(sum = c("spatialized", "total", "diff"), rbind(sum.p.1A4ai, total.1A4ai, data.frame(sum.p.1A4ai == total.1A4ai)-1)) %>%
  datatable(., caption = 'Table 3: Summary differences after spatialisation',
            options = list(pageLength = 5)
  )

# Export data to geopackage file
# st_write(p.1A4ai, dsn="Products/1A4 - Residential-Tertiary/1A4ai.gpkg", layer='1A4ai')